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TITLE OF THE INVENTION 
APPARATUS MSID IflETHOD FOR PERFORMING TIME DELAY ESTIMATION OF 
SIGNALS PROPAGATING THROUGH AN ENVIRONMENT 

10 

CROSS REFERENCE TO RELATED APPLICATIONS 
This application claims priority of U.S. Provisional Patent 
Application No. 60/494,358 filed August 12, 2003 entitled ECHO 
DELAY ESTIMATES FROM MULTIPLE SONAR PINGS. 

15 

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR 

DEVELOPMENT 

This invention was made with government support under U.S. 
Government Contract Nos. ARO DAAD 19-02-1-0403 and ONR N00012-02- 
20 C-02960. The government has certain rights in the invention. 

BACKGROUND OF THE INVENTION 
The present application relates generally to signal 
processing, and more specifically to' systems and methods of 
25 increasing the accuracy of time delay estimates of signals 
propagating through an environment. 

Various industrial and scientific techniques require 
accurate estimations of time delays of signals propagating through 
an environment such as an underwater environment, a soil 
30 environment, or an environment comprising living tissue. For 
example, in an underwater environment^ sonar systems may be 
employed to estimate time delays of sonar pulses reflected from an 
object or target to estimate a distance to the target (also known 
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as estimating the range of the target) . Conventional systems for 
performing sonar range estimation typically transmit one or more, 
sonar pulses (^^pings") comprising sonic or supersonic pressure 
waves toward a selected target, and receive one or more sonar 
5 pulses reflected from the target. Such reflected sonar pulses 
(^^ echoes" or returns'') may include a significant amount of 
background noise and/or other interfering signals in addition to 
the reflected sonar signals of interest. For example, a 
conventional sonar system may comprise a coherent receiver 

10 including a cross correlator configured to receive the echo and a 
representation of the transmitted sonar pulse or ping, which are 
cross-correlated within the coherent receiver to generate a peak 
cross correlation value. The conventional sonar system typically 
compares the peak cross correlation value to a predeteimined 

15 threshold value. If the cross correlation value is greater than 
the predetermined threshold value, then the reflected sonar signal 
of interest has been successfully detected. The conventional 
sonar system may then utilize the cross correlation peak to obtain 
a measure of the range of the target. 

20 One drawback of the above-described conventional sonar 

system is that the level of background noise and/or other 
interfering signals contained within the echo or return may be 
sufficient to cause the reflected sonar signal to go undetected or 
to be falsely detected, thereby causing the cross correlator to 

25 produce inaccurate range measurements. Such inaccurate range 
measurements are likely to occur in low signal-to-noise ratio 
(SNR) sonar environments, in which the noise power within the echo 
may be comparable to or greater than the reflected signal power. 
This can be problematic in sonar ranging systems because a 

30 reduction in the measurement accuracy of the cross correlator 
typically leads to a concomitant reduction in sonar ranging 
accuracy. 
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Prior attempts to increase the accuracy of sonar ranging 
have included filtering out at least some of the background noise 
before providing the echoes to the cross correlator. However, 
such attempts have generally not workeci well enough to allow 
successful detection of reflected sonaar signals and accurate 
estimation of range in low SNR sonar environments . This is due, • 
at least in part, to the fact that sonar systems typically receive 
echoes that include various types of noise from a variety of 
different noise sources. For example, a sonar system may transmit 
pings through a medium such as water from a ship or submarine that 
produces noise across a wide frequency range. Further, other 
ships, submarines, or structures producing noise across wide 
frequency ranges may be within the vicini-ty of the sonar system. 
Moreover, the natural interaction of the water and objects within 
the water including the selected target m.ay produce a substantial 
amount of ambient noise. 

In addition, sonar ranging systems may receive echoes from 
multiple selected (and unselected) targets , each target having its 
own associated noise level, and it may be desirable to determine 
the noise level and range of each target separately. Such noise 
associated with multiple targets may be stationary or non- 
stationary, linear or nonlinear, or additive or non-additive. 
Further, at least some of the backgrouncd noise may result from 
reverberations and/or random signal distorrtions of the ping and/or 
echo, and therefore the noise level ancd its structure may be 
significantly affected by the transmitted sonar signal. However, 
conventional sonar systems are generally incapable of accurately 
estimating noise levels and target ranges in the presence of non- 
stationary, nonlinear, non-additive, and/or signal-dependent 
noise. 

Moreover, the density and temperature of the transmission 
medium (e.g., water) and the frequency of the 
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transmitted/reflected sonar signals may affect the decay rate of 
the sonar pulse propagating through the medium. In addition, the 
absorption of certain frequencies of the ping by the target may 
affect the strength of the resulting echo. However, conventional 
5 sonar systems are generally incapable of fully compensating for 
such factors when called upon to generate accurate noise and range 
estimates . 

It would therefore be desirable to have a system and method 
of increasing the accuracy of time delay estimates of signals 
10 propagating through an environment that avoids the drawbacks of 
the above-described systems and methods. Such a system would have 
increased resilience to noise, thereby allowing an increase in the 
operating range of the system and/or a decrease in the power level 
of signals employed by the system. 

15 

BRIEF SUMMARY OF THE INVENTION 
In accordance with the present invention, a system and 
method are provided for increasing the accuracy of time delay 
estimates of signals propagating through an environment. The 

20 presently disclosed system and method achieve such increased 
accuracy in time delay estimation Joy employing multiple 
transmitted signals and/or multiple received signals. Further, in 
the event a degree of noise accompanies the received signals, at 
least some of the noise is non-correlated. 

25 In one embodiment, the system includes one or more sensors 

for receiving a plurality of signals, and a time delay estimator 
for measuring time delays between multiple pairs of the plurality 
of signals, thereby generating time delay estimation data from the 
measured time delays. At least some of the multiple pairs of 

30 signals are received and measured at different points in time. 
The system further includes a data analyzer for analyzing the time 
delay estimation data, for generating a statistical distribution 
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of the time delay estimates from the tint^ delay estimation data^ 
and for calculating at least one of the m-^an^ the median^ and the 
mode of the time delay estimation distribua.tion. By increasing the 
number of signals employed by the system^ the accuracy of the time 

5 delay estimation is increased. Further^ IzDy calculating the median 
or the mode of the distribution of time <lelay estimations^ noise 
tolerance is improved. The plurality of signals may comprise one 
of sonar signals, seismic signals, ultrsisonic signals, acoustic 
signals, electromagnetic signals, or any , other suitable type of 

.0 signals. 

In another embodiment, the system. and method employ an 
iterative process in which the distir -xbution of time delay 
estimates is determined, and an initial statistical estimate of 
time delay is calculated from the distribuition. Next, a first set 

.5 of boundaries of the time delay estJ_ination distribution is 
determined, and time delay estimates dispc:>.sed outside of the first 
set of boundaries (^^ outliers" ) are remove <d from the distribution. 
The statistical estimate of time delay ±-s then recalculated and 
compared to the initial statistical estixmate. In the event the 

:0 difference between the initial statists ical estimate and the 
recalculated statistical estimate is greater than or equal to a 
predetermined threshold value, a nextzz set of distribution 
boundaries is determined. Next, one or amore additional outliers 
are removed from the distribution, and thst statistical estimate of 

:5 time delay is recalculated and compared to the prior statistical 
estimate. In the event the difference tz>etween the recalculated 
statistical estimate and the prior static stical estimate is less 
than the predetermined threshold value, the final statistical 
estimate is used to increase the acr:curacy of time delay 

JO estimation. In alternative embodiments, -the statistical estimate 
of time delay comprises the mean, the mecaian, or the mode of the 
distribution of time delay estimations. 
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In general^ the mean of a distribution of multiple 
observations is obtained by multiplying the value of each 
observation by its probability and summing the resulting products . 
In general, the median of a distribution of multiple observations 
5 is obtained by ranking the values of the observations in order 
from smallest to largest and taking the central value. In 
general, the mode of a distribution of multiple observations is 
obtained by taking the most frequent value of the observations. 

Other features, functions, and aspects of the invention will 
10 be evident from the Detailed Description of the Invention that 
follows • 

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS 
The invention will be more fully understood with reference 
15 to the following Detailed Description of the Invention in 
conjunction with the drawings of which: 

Fig. la is a block diagram of a system for performing time 
delay estimation of signals propagating through an environment 
according to the present invention; 
20 Fig. lb is a block diagram of an alternative embodiment of 

the system of Fig. la; 

Figs. 2a-2b are diagrams of ambiguity functions illustrating 
the effect of noise level on the variability of cross correlation 
peaks; 

25 Fig. 3a is a diagram illustrating peak variability as a 

function of signal-to-noise ratio and center frequency for a 

plurality of frequency bands; 

Fig. 3b is a diagram illustrating a performance curve 

derived from the diagram of Fig. 3b; 
30 Fig. 4 is a flow diagram illustrating a method of operating 

a first embodiment of the system of Fig. La; 
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Fig. 5a is a diagram illustrating characteristics of one of 
multiple pings transmittable by the system of Fig. la; 

Fig. 5b is a diagram illustrating characteristics of 
successive echoes resulting from the ping of Fig. 5a; 
5 Figs. 6a-6c are diagrams illustrating the successive echoes 

of Fig. 5b after being preprocessed by a signal processor included 
in the system of Fig. la; 

Figs. 7a-7c are diagrams illustrating probability density 
functions for pulse arrival time delays (PATDs) corresponding to 
10 the preprocessed echoes of Figs. 6a-6c; 

Fig. 8 is a diagram illustrating a distribution of PATDs for 
successive echoes returning from a target cylinder filled with 
kerosene; 

Figs. 9a-9c are diagrams illustrating probability density 
15 functions for mean pulse arrival time delays (MPATDs) 
corresponding to the preprocessed echoes of Figs. 6a-6c; 

Fig. 10a is a conceptual model of a noiseless cross 
correlation function corresponding to a single echo/ping pair; 

Fig. * '^lOb is a diagram illustrating the probability of 
20 selecting a given bin in the cross correlation function of Fig. 
lOa; 

Fig. 11a is a second model of a cross correlation function 
corresponding to a single echo/ping pair; 

Fig. lib is a diagram illustrating the probability of 
25 selecting a given bin in the cross correlation function of Fig. 
11a; 

Fig. 12 is a third model of a cross correlation function 
corresponding to multiple echo /ping pairs; 

Fig. 13a-13d are diagrams of simulated distributions of echo 
30 delay estimates for different noise levels in the environment of 
Fig. la; 
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Figs. 14a-14ci are diagrams illustrating peak variability as 
a function of signal-to-noise ratio for different numbers of 
pings; and 

Fig. 15 is a diagram of performance curves comprising 
5 composites of the peak variability curves of Figs. 14a-14d, 
showing accuracy breakpoints versus the number of pings. 

DETAILED DESCRIPTION OF THE INVENTION 
U.S. Provisional Patent Application No. 60/494^358 filed 

10 August 12 r 2003 entitled ECHO DELAY ESTIMATES FROM MULTIPLE SONAR 
PINGS is incorporated herein by reference. 

Systems and methods are disclosed for increasing the- 
accuracy of time delay estimates of signals propagating through an 
environment. The presently disclosed system and method increase 

15 time delay estimation accuracy by using time delay estimates 
generated from multiple signals including varying levels of non- 
correlated noise to form a statistical distribution of the time 
delay estimates. In one embodiment^ the system and method employ 
an iterative process in which multiple signals are preprocessed, a 

20 distribution of time delay estimates of the signals is determined/: 
and the mean of the distribution is calculated. In alternative 
embodiments, the median or the mode of the distribution of time 
delay estimates is calculated. By calculating the mean^ the 
median, or the mode of a distribution of time delay estimates, 

25 information from multiple signals may be employed to increase the. 
accuracy of the time delay estimation. In addition, by 

calculating the median or the mode of the distribution of time 
delay estimates, the accuracy of time delay estimation is 
increased while improving the noise tolerance of the system. 

30 Fig. la depicts an illustrative embodiment of a system 100 

for estimating time delays of signals propagating through an 
environment, in accordance with the present invention. In the 
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illustrated embodiment, the system 100 includes one . or more 
sensors 102. 1-102. n, a transmitter 103, a receiver 104, a time 
delay estimator 108 including a signal processor 106, and a data 
analyzer 110. It is noted that the illustrative embodiment of the 
5 system 100 described herein is suitable for performing time delay 
estimation in (1) a fluid environment, e.g., an air environment, 
or an underwater environment for marine exploration, (2) an earth 
environment for seismic exploration, (3) an environment comprising 
living tissue for medical ultrasound, or any other suitable 

10 environment. The signals propagating through the various 

environments may therefore comprise sonar signals, seismic 
signals, ultrasonic signals, acoustic signals, electromagnetic 
signals, or any other suitable type of signals. It should also be 
understood that the presently disclosed system 100 may be adapted 

15 for use in radar systems, microwave systems, laser systems, or any 
other suitable system. 

In the presently disclosed embodiment, the transmitter 103 
is configured to transmit one or more sonar pulses (^^pings'') 
through a transmission mediiom such as water. The pings travel 

20 through the water until they strike an object or target 112 in the 
water, which returns one or more reflected sonar pulses (^^ echoes'' 
or ^^returns") toward the sonar sensors 102. 1-102. n. For example, 
the sonar sensors 102. 1-102. n may comprise one or more hydrophone 
sensors. Each one of the sensors 102. 1-102. n is configured to 

25 sense the echoes, and to provide signals representative of the 
echoes to the sonar receiver 104. In turn, the receiver 104 
provides indications of the echoes to the time delay estimator 
108. 

In the illustrated embodiment, the time delay estimator 108. 
30 receives the echo indications from the receiver 104, and receives 
representations of the pings transmitted by the sonar transmitter 
103. For example, the time delay estimator 108 may comprise a 
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cross correlator^ or any other suitable device or technique for 
estimating the time delay of signals. In the presently disclosed 
embodiment^ the time delay estimator 108 comprises a cross 
correlator^ which is configured to perform multiple cross 
5 correlation operations on the echoes and pings. Specifically/ the 
time delay estimator 108 cross-correlates each echo/ping pair, and 
provides cross correlation output data to the data analyzer 110 , 
which is operative to analyze the data to make multiple echo delay 
estimates by determining the variability of cross correlation 

10 peaks, and to generate a distribution of the echo delay estimates. 
The preprocessing of the echoes and the generation of the echo 
delay estimation distributions are described below. It should be 
understood that the system 100 may comprise an active system 
capable of estimating time delays of signals using multiple 

15 transmitted signal/return signal pairs , or a passive system 
capable of estimating time delays using successive received 
signals. In the event the system 100 comprises a passive system, 
the transmitter 103 may be removed from the system. 

The operation of the presently disclosed sonar system 100 

20 will be better understood by reference to the following analysis.:. 
In general, the cross correlation of an echo and a ping may be 
expressed as 

\|/e O \|/p(x) = J\|/e(t)\]/p(t+T) dt => 
25 = J\|/p(t)\l/p(t+T+To)dt + J\|/p(t)T|(t+T)dt, (1) 

in which the first term ^^J\|/p (t ) \|/p (t+x+Xo) dt" is the auto- 
correlation of the ping centered at time Xo (e.g., Xo=0) , the 
second term ''J\)/p (t ) t] (t+x) dt" is representative of band-limited 
30 white noise with frequency limits defined by the spectrum of the 
ping, and the integration operation in each term is performed from 
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-00 to +00. Figs. 2a-2b depict representative ambiguity functions 
201-203 that may be employed to describe the output provided by 
the time delay estimator 108 (see Fig. la) comprising a cross 
correlator- Because the cross correlator operates on pairs of 
5 echoes and pings, it is understood that an ambiguity function may 
be constructed corresponding to each echo/ping pair. 

As shown in Figs. 2a-2b;. the ambiguity functions 201-203 are 
expressed as functions of sonar pulse amplitude (vertical axis, 
dB) and delay time (horizontal axis, seconds), which is 

10 proportional to sonar range in a sonar ranging system. . 
Specifically, the ambiguity functions 201-203 correspond to the 
cross correlation of respective echo/ping pairs having 
approximately the same frequency range but different center 
frequency fc (i.e., mean integrated frequency). For example, the 

15 ambiguity function 201 corresponds to the cross correlation of a 
first echo/ping pair having a low center frequency fcl, the 
ambiguity function 202 corresponds to the cross correlation of a 
second echo/ping pair having an intermediate center frequency fc2, 
and the ambiguity function 203 corresponds to the cross 

20 correlation of a third echo/ping pair having a high center 
frequency fc3. Fig. 2a depicts a detailed view of the main lobes 
of the ambiguity functions 201-203, and Fig. 2b depicts the main 
lobes and side lobes of the ambiguity functions 201-203. Each one 
of the ambiguity functions 201-203 comprises a respective peak 

25 value, which is indicative of the range of the target returning 
the echo in a sonar ranging system. 

In high SNR sonar environments (i.e., when the noise level 
is low) , the peak of the ambiguity function is generally located 
at the main lobe of the function. In this case, the peaks of the 

30 ambiguity functions 201-203 are regarded as having low ambiguity, 
and may be located within the width of the main lobes of the 
functions at about time Xof as illustrated by the vertical line of 
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Pig. 2a. It is appreciated that the time tq corresponds to the 
actual range of the target. The effect of the low level of noise 
in the sonar environment is to jitter the position of the peak 
around the time To- To a first approximation^ the magnitude of 
5 this jitter (also known as peak variability) is relatively low, 
e.g., the peak variability is typically less than the width of the 
main lobes, as illustrated by the horizontal lines of Fig. 2a. 
The lengths of the horizontal lines of Fig. 2a are indicative of 
the levels of peak variability associated with the respective 
10 ambiguity functions 201-203. In the illustrated embodiment, the 
lowest peak variability is associated with the ambiguity function 
203 (high center frequency fc3) , and the highest peak variability 
is associated with the ambiguity function 201 (low center 
frequency fcl) . 

15 In low SNR sonar environments (i.e., when the noise level is 

high, for example, when the noise level is of the order of the 
difference between the amplitudes of the main lobe and the first 
side lobe) , the peak of the ambiguity function may not be located 
within the main lobe of the function, but instead may be located 

20 at one of the side lobes. In this case, the peaks of the 
ambiguity functions 201-203 are regarded as having high ambiguity, 
and may be located (1) within the width of a side lobe at about 
time X-2 for the function 203, (2) within the width of a side lobe 
at about time for the function 202, and (3) within the width 

25 of a side lobe at about time t-3 for the function . 201, as 
illustrated by the vertical lines 203a, 202a, and 201a, 
respectively, of Fig. 2b. The effect of the high level of noise* 
in the sonar environment is to significantly increase the peak 
variability, thereby increasing the potential error in sonar 

30 ranging. The horizontal lines in Fig. 2b illustrate the potential 
error in sonar ranging that can result from such high noise 
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levels. The peak variability associated with the ambiguity 

functions 201-203 is indicative of the accuracy of echo delay 

estimates for the respective echo/ping pairs. 

Fig. 3a depicts peak variability as a function of SNR (dB) 
5 and center frequency fc for a plurality of frequency bands. As 

shown in Fig, 3a;. peak variability is expressed in terms of root 

mean square error (RMSE, seconds) , which is a temporal 

representation of the potential error in echo delay estimation. 

Further,, the center frequencies fc of the frequency bands are 
10 equal to 12 kHz, 25 kHz, 37 kHz, 50 kHz, 62 kHz, 75 kHz, and 87 

kHz, respectively, and the centralized root mean square bandwidth 

BcBMs of the pings is fixed at 2 . 1 kHz . 

For example, a plurality of peak variability curves 301-307 

(see Fig- 3a) may be obtained via Monte Carlo simulations. 
15 Specifically, sonar pings may be expressed as cosine packets of 

the form 

M/cT,T,(t) - Ka,i,exp(-t^/2{STD)2)cos(27CTit), (2) 

20 in which '''ti''' is the center frequency, ^^STD" is the standard 
deviation of a peak location in time which controls the spread in 
time of the ping and its frequency bandwidth, and ^"Ka^i^" is a 
normalization factor such that 

25 Va,^(t)dt - 1, (3) 

in which the integration operation is performed from -oo to 4-oo. 
Further, white noise may be added to the pings to generate noisy 
echoes for simulation purposes, and a temporal indication of the 
30 echo delay estimate may be computed as the time corresponding to 
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the maximiam amplitude of the cross correlation between the echo 
and the ping. 

As shown in Fig, 3a, each one of the simulation curves 301- 
307 is approximately linear within a first SNR range of about 35- 
5 50 dB (see also region I of Fig, 3b) . Further, for each curve 
301-307;- there is a sharp transition from lower RMSE levels to 
higher RMSE levels within a second SNR • range of about . 15-35 dB 
(see also region II of Fig. 3b), thereby indicating significant 
increases in peak variability. Within a third SNR range of about' 

10 5-15 dB (see also region III of Fig. 3b), the curves 301-307 are 
again approximately linear. It is noted that the curve 308 
depicted in Fig. 3b is a performance curve comprising a composite 
of the peak variability curves 301-307, including breakpoints 1-9 
(see Fig. 3a). Accordingly^ as the SNR decreases (i.e., as the 

15 noise level increases) , the RMSE levels gradually increase within 
region I until sharp transitions occur from lower RMSE levels to 
significantly higher RMSE levels within region II - the RMSE 
levels then continue to increase more rapidly within region III. 
It is noted that within region III, the sonar range resolution 

20 falls sharply until the sonar is ineffective and the target is 
considered to be out-of-range. 

Specifically, within region I, the simulation curves 301-307- 
approximately track a line 310 (see Fig. 3b), which may be defined 
as 

25 

STD = (27cBRMsd)'"S (4) 

in which ^^STD" is the standard deviation of a peak location in 
time and is proportional to the RMSE, ^^BpMs" is the root mean 
30 square bandwidth of the ping, and is the SNR. A derivation of 

equation (4) is described in ProJbaJbiiitY and Information Theory 
with Applications to Radar, P.M. Woodward, New York, McGraw-Hill 
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5 



Book Company, Inc., copyright 1953, which is incorporated herein 
by reference. It is noted that Brms niay be expressed as 

Bkms = (Jf^PsD(f)df)^/^ (5) 

in which ^'PsD(f)'' is the power spectral density of the ping, and 
the integration operation is performed from 0 to +oo. Further, d 
may be expressed as 

10 d = (2E/No)^^S (6) 

in which ^^E" is the total energy of the echo, and ^^Nq" is the. 
spectral density of the noise. Accordingly, 

15 SNR(dB) = 201ogiod. (7) 

Moreover, following the sharp transitions from lower RMSE 
levels to higher RMSE levels within region II (see Fig. 3b), the 
simulation curves 301-307 approximately track a line 312 (see Fig. 
20 3b) , which may be defined as 

STD = (27CBcRMsd)"^, (8) 

in which ^^Bcrms" is the centralized root mean square bandwidth of 
25 the transmitted pulse. The RMSE levels continue to increase at a 
faster rate within region III (see Fig. 3b) . It is noted that 
BcKMs ni^y be expressed as 



BcBMs - (J(f-fc)2psD(f)df)^/2^ (9) 



30 
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in which ^^fc" is the center frequency of the ping, and the 
integration operation is performed from 0 to +00. it is further 
noted that fc may be expressed as 

5 f c = JfPsD(f)df, (10) 

in which the integration operation is performed from 0 to +00. 
Moreover, the root mean square bandwidth may be expressed as 

10 BrmS^ = BcRMS^ + fc^. (11) 

Accordingly, in the event the center frequency fc is much larger 
than the centralized root mean square bandwidth Bcrms/ 

15 Brms « fc. (12) 

The behavior of the simulation curves 301-307 within region 
I (see Figs. 3a-3b) is characteristic of the performance of a 
coherent"' receiver, which estimates the echo delay relative to a 

20 peak of the ambiguity function within the width of the function's 
main lobe. The behavior of the curves 301-307 after their sharp 
transitions from lower RMSE levels to higher RMSE levels within 
region II (see Figs. 3a-3b) is characteristic of the perfoimance 
of a '^semi-coherent" receiver, which estimates echo delay relative 

25 to the peak of the envelope of the ambiguity function. As 
illustrated in Fig. 3a, echo delay estimates provided by the semi- 
coherent receiver have associated errors (RMSE) that are 
significantly higher than the errors associated with the echo 
delay estimates of the coherent receiver. 

30 In a first embodiment of the sonar system 100, the accuracy 

of echo delay estimation is increased via a method employing 
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multiple echo/pulse pairs in an active sonar system, or successive 
echoes in a passive sonar system. The method includes estimating 
the distribution of echo delays corresponding to multiple pings, 
and eliminating echoes that provide echo delay estimations that 
5 lie outside defined boundaries of the distribution. The remaining 
echoes are then employed to improve the echo delay estimation. 
The method of operation of this first embodiment is described- 
below by reference to Figs. 1 and 4. 

As depicted in step 402, multiple echoes are received by the 

10 sensors 102. 1-102. n, and are provided in turn to the receiver 104 
and to the time delay estimator 108, which includes the signal 
processor 106. The multiple echoes are then preprocessed, as 
depicted in step 404, by the signal processor 106. In the event 
the time delay estimator 108 comprises a cross correlator, 

15 multiple echo/ping pairs are cross correlated, as depicted in step' 
406, and cross correlation output data is provided to the data 
analyzer 110. It is noted that at least some of the respective 
echo/ping pairs are received and cross correlated at different 
points in time. Respective echo delay estimations are then 

20 determined, as depicted in step 408, and a distribution of the 
echo delay estimates is generated, as depicted in step 410, by the 
data analyzer 110. Next, an initial mean of the distribution of 
echo delay estimates is calculated, as depicted in step 412, by 
the data analyzer 110. A first set of boundaries of the 

25 distribution is then determined, as depicted in step 414, and echo 
delay estimates lying outside the first set of boundaries 
('^outliers'') are removed, as depicted in step 416, by the data 
analyzer 110. Next, the mean of the distribution is recalculated, 
as depicted in step 418, and the recalculated mean is compared to 

30 the initial mean, as depicted in step 420, by the data analyzer 
110. In the event the difference between the recalculated mean 
and the initial mean is greater than or equal to a predetermined 
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threshold value, as depicted in step 422, a next set of 
distribution boundaries is determined, as depicted in step 414, by 
the data analyzer 110. Next, additional outliers are removed from 
the distribution based on the next set of distribution boundaries, 
5 as depicted in step 416, by the data analyzer 110. The mean of 
the distribution is then recalculated, as depicted in step 418, 
and compared to the prior calculated mean, as depicted in step 
420, by the data analyzer 110. In the event the difference 
between the recalculated mean and the prior mean is less than the 

10 predetermined threshold value, as depicted in step 422, the 
process ends and the final recalculated mean of the echo delay- 
estimation distribution is used to estimate the echo delay with 
increased accuracy. In the event the difference between the 
recalculated mean and the prior mean is not less than the 

15 predetermined threshold value, the process loops back to step 414 
to determine a new set of distribution boundaries . 

The method of Fig. 4 will be better understood by reference 
to the following illustrative example. In this example, increased 
accuracy in echo delay estimation is illustrated via improved. 

20 discrimination of a target comprising a cylindrical container 
filled with different liquids such as fresh water, saline water, 
glycerol, or kerosene- Multiple pings are transmitted toward the 
target, and returning echoes are recorded by one or more 
hydrophone sensors. Because physical properties such as density. 

25 and compressibility of the liquids differ from each other, a 
liquid inside the cylindrical target can be identified based on 
the sound velocity within the liquid. The multiple pings 
transmitted toward the target are broadband and have a prominent 
peak 502 in the time domain, as depicted in Fig. 5a. When each 

30 ping penetrates the cylindrical target, the peak is maintained 
after being bounced from the front and the back of the cylinder. 
Two peaks 503-504 can therefore be detected in the successive 
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echoes, as depicted in Fig. 5b, and the time difference between 
the two peaks 503-504 can be measured. Accordingly, 
discrimination of the liquid within the target can be achieved by 
analyzing the temporal differences between successive echoes, 
which are referred to herein as the peak arrival time differences 
(PATDs) . 

Specifically, the PATD is defined as the difference between 
the peak arrival times of successive echoes. For example, the 
peak arrival times of successive echoes may be expressed as 

tmax.a = arg max y(t) (13) 
tmax^b = arg max y' (t) , (14) 

in which ''y(t)" and ^"y' (t)'' designate preprocessed sonar echoes. 
It is noted that y' (t) is generated with the peak near t^ax^a 
removed. Accordingly, the PATD may be expressed as 

PATD s |t„^x,a - troax.bl. (15) 

It is noted that sonar pings may be subject to phase 
inversion when they are reflected by low density media such as 
kerosene. Further, the detection of peaks in echoes in the time 
domain is often sensitive to noise. Moreover, the wave shape of 
the echoes may change as they propagate in time due to dispersion. 
As a result, there may be a significant amount of uncertainty in 
the temporal locations of peaks in the echoes. 

For this reason, the echoes are preprocessed to make the 
deteirmination of the PATD more stable. In the presently disclosed" 
embodiment, the signal processor 106 (see Fig. la) employs a 
plurality of preprocessing techniques. In a first preprocessing 
technique, the absolute values of the echoes are taken to prevent 
peak loss by phase inversion, i.e., 
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25 



ya(t) = |x(t) I . (16) 

In a second preprocessing technique;- a matched filter such 
as a linear filter is employed in the receiver 104 (see Fig. la) . 
to enhance signal detectability by increasing the signal-to-noise 
ratio (SNR) . In the preferred embodiment, the time-reversed 
pinging signal yields the maximum SNR, and the resulting output 
signal is expressed as 

ym(t) = f^x(t')s(t'-t)dt'^xit)*s(t% (17) 



in which ^^s(t)'^ and ^^x(t)" represent the pings and the echoes, 
respectively . 

15 In a third preprocessing technique, the wave envelope peak 

is detected to reduce temporal uncertainty. In the preferred 
embodiment, instantaneous envelope detection (lED) is employed to 
detect abrupt changes in the wave amplitude. The instantaneous 
envelope is the amplitude of an analytic signal, of which the real 

20 part is the recorded signal and the imaginary part is a Hilbert 
transform of the signal. For example, the Hilbert transform x(t) 
of a signal x(t) may be expressed as 

ji(0 = x(t)*l- = ^f^^dt\ (18) 
nt jt '-^ t-t' 

Accordingly, the envelope function ye(t) may be expressed as 

y^it)=\x(t)+ix(t)\ (19) 
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Figs. 6a-6c depict the successive echoes 503-504 of Fig. 5b, 
and preprocessed versions 603a"-603Cr 604a-604c of the successive 
echoes 503-504. Fig. 6a depicts the echoes 503-504 preprocessed 
to obtain the absolute values 603a-604a of the echoes. Fig. 6b 
5 depicts the echoes 503-504 preprocessed by a matched filter to 
obtain the filtered versions 603b- 604b of the echoes, and Fig. 6c 
depicts the echoes 503-504 preprocessed to obtain the Hilbert 
transforms 603c-604c of the echoes. It is understood that the 
time delay estimator 108 (see Fig. la) may determine the PATD' s 
10 directly from the preprocessed signals 603a-604a, 603b-604b, and 
603c-604c. 

Figs. 7a-7c depict probability density functions (PDFs) of 
the PATDs calculated after preprocessing the echoes by the 
absolute value technique (see Fig. 7a), the matched filtering 

15 technique (see Fig. 7b) , and the instantaneous envelope detection 
technique (see Fig. 7c) . As shown in Figs. 7a-7c, the PDFs 702a- 
702c correspond to the target cylinder filled with fresh water, 
the PDFs 704a-704c correspond to the target cylinder filled with 
saline water, the PDFs 706a-706c correspond to the target cylinder 

20 filled with glycerol, and the PDFs 708a-708c correspond to the 
target cylinder filled with kerosene. As shown in Fig.- 7a, the 
PDF 708a corresponding to kerosene has two peaks next to each 
other. This is because the absolute value preprocessing technique 
discards the phase information, and the PDF peaks therefore 

25 include both positive and negative peaks. As shown in Fig. 7b, 
the PDF 704b corresponding to saline includes multiple peaks. 
Specifically, the saline peaks indicate periodic error, which may 
be attributed to the coherence of sound waves between the first 
and second peaks. Such coherence produces a resonant standing 

30 wave inside the target cylinder filled with saline water. Because 
the matched filter tends to reinforce the carrier frequency, the 
resonance is emphasized. As shown in Fig. 7c, the PDFs 702c, 
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704c, 706c, and 708c detennined after preprocessing with the 
instantaneous envelope detection technique are relatively stable 
and robust. 

Fig. 8 depicts an exemplary distribution 802 of PATDs for 
5 successive echoes returning from the target cylinder filled with 
kerosene • As shown in Fig. 8, a set of boundaries of the PATD 
distribution is indicated by reference numerals 804 and 806. In 
this example, the largest PATD distribution peaks of four classes 
are included within the boundaries 804 and 806, and the PATD 

10 distribution peaks lying outside of the boundaries (''outliers") 
are removed. Such outliers may be similarly removed from the PATD 
distributions corresponding to the target cylinder filled with 
fresh water, saline water, and glycerol. Next, the mean values of- 
the PATD distributions are calculated without the outliers to 

15 obtain the mean PATD distributions (MPATDs) for fresh water, 
saline water, and glycerol, and kerosene. 

Figs. 9a-9c depict the probability density functions (PDFs) 
of the MPATDs calculated after preprocessing by the absolute value 
technique (see Fig. 9a), the matched filtering technique (see Fig. 

20 9b), and the instantaneous envelope detection technique (see Fig. 
9c). As shown in Figs. 9a-9c, the PDFs 902a-902c correspond to 
the target cylinder filled with fresh water, the PDFs 904a-904c 
correspond to the target cylinder filled with saline water, the 
PDFs 906a-906c correspond to the target cylinder filled with 

25 glycerol, and the PDFs 908a-908c correspond to the target cylinder 
filled with kerosene. As indicated by the narrower widths of the 
peaks in the MPATD PDFs for fresh water, saline water, glycerol, 
and kerosene, the discrimination of each liquid is improved. It 
is understood that a new set of boundaries of the MPATD 

30 distribution may be determined to remove additional outliers, and 
the mean values of the distributions may be recalculated to 
further improve target discrimination. 
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In a second embodiment of the sonar system 100 (see Fig. 
la) , the accuracy of echo delay estimates is increased via a 
method that includes calculating the mode of echo delay estimation 
distributions generated from raw or preprocessed sonar signals. 
As described below^ the method of this second embodiment not only 
increases the accuracy of echo delay estimates ,r but it also 
increases the noise tolerance of the system 100. 

The method of the second embodiment of the sonar system 100 
will be better understood by reference to the following analysis. 
Fig. 10a depicts a conceptual model of a noiseless cross 
correlation function 1002 corresponding to a single echo/ping 
pair. As shown in Fig. lOa^ the function 1002 is a delta (5) 
function that is zero everywhere within an a priori window of 
length L, except at the time location (^^bin'') t=0 when it is equal 
to the energy of the ping. It is noted that this conceptual model 
corresponds to a sonar ping having infinite bandwidth. When white 
noise is added to the echo, the cross correlation function 1002 
has a Gaussian distribution with multidimensional centers at zero 
for all values that are outside of the bin t— 0, and a center at 
a=E>0 for the value of the function 1002 within the bin t=0, in 
which E is the energy of the echo. If xi is the value of the 
function 1002 within each bin t=i, and xo is the value within the 
bin t=0, then a probability density function (PDF) may be 
expressed as 

P(Xr>X.>-X.) = -^^^e . (20) 
v=[a,0,0,...,0]^ 

Using the Gaussian distribution, the probability "a" of the event 
xl>xi (i^fel) , which corresponds to the "correct" echo delay 
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estimation for a given noise level a (l.e.r for a given SNR) ^ may 
be expressed as 

00 X, X. xt ^ i|x-v[|' 

P(X^>X^,i:stl) = a=^ i dx, I dx:, i dx^.,. j dx^ rL= — e ^-^^ => 

-00 -00 -<o -00 (<Tv2;r)^ 

a= J- 1 e~^"^^\l+effix)]''-'dx. (21) 



It is noted that the probability a expressed in equation (21) 
depends on " a which is proportional to the signal-to-noise 

10 ratio (SNR) • A probability of error "^^p" may therefore be 
expressed as 1-a^ which represents the probability that the 
amplitude of at. least one peak disposed outside of the correct bin 
is larger than the amplitude of the peak disposed within the 
correct bin- Fig. 10b depicts the probability a of selecting a 

15 given bin in the cross correlation function 1002. 

Fig. 11a depicts a more practical model in which a cross 
correlation function 1102 of an echo/ping pair comprises a uniform 
distribution within an interval ^^A" . It is noted that this second 
model corresponds to a sonar ping having a finite bandwidth. In 

20 the model of Fig. 11a, the cross correlation function 1102 is 
approximated by a piecewise constant function having an amplitude 
^^a'' within the central interval Ia and zero elsewhere. For 
example / the length L of the a priori window of the cross 
correlation function 1102 may correspond to sonar range. In the 

25 event the a priori window has a length of 2L and a sampling 
frequency of ^^fs"' is employed, the total number of echo delay 
estimations may be expressed as 
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N = Na + No = 2L • fs, (22) 

in which ''^Na=A«fs"' estimates are disposed within the central or 
correct" bin t==0 and ''^No"' estimates are disposed outside of the 
5 correct bin but within the a priori window. A random vector may 
be defined in which the first Na random variables correspond to 
the amplitudes of the estimates within the correct bin, and the 
last No random variables correspond to the amplitudes of the 
estimates outside of the correct bin, in which Na and No are 
10 integers. As indicated above, for white noise, the joint 
probability density function for the vector of random variables, 
may be expressed as 



ll^-g|i^ 



15 



20 



(crV^^r) (23) 



V = [a, a,..., a, 0,0,...,0] 



T 



0 



The probability of time delay estimation within the correct bin of 
the cross correlation function 1102 may therefore be expressed as 

(24) 

2 -ylTT"^ /V2cr 



Accordingly, for a given noise level ^^o" (i.e., for a given SNR) 
the probability a that the correct bin is selected may be 
expressed as 
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<^ = tJt^ ^'\i^eff{x--^)f^^'[l + e^^^ (25) 

Fig. lib depicts the probability a of selecting a given bin 
in the cross correlation function 1102. 

As described above with reference to Fig. 3a, the peak 
variability curves 301-307 of Fig. 3a include sharp breakpoints 
between the curve regions corresponding to the coherent, receiver 
and the curve regions corresponding to the semi -coherent receiver. 
As explained above, these breakpoints in the curves 301-307 
indicate that echo delay estimation accuracy is sharply reduced as 
the SNR of the environment decreases. The accuracy breakpoints in 
peak variability curves are further analyzed below. 

In this analysis, a random vector whose probability 
distribution is given by equation (24) is designated as ^^T^' . For 
an SNR value that is greater than a given value SNRo, the 
probability of an echo delay estimate falling within the correct 
bin is greater than a given probability ao. The standard 
deviation STD of the distribution of echo delay estimation within 
the correct bin is designated as ^^STDa", and the standard 
deviation STD of the distribution outside of the correct bin is 
designated as ^^STDo''. The distribution T whose cumulative 
function is expressed by equation (25) is sampled. For n" 
estimations, a fraction of the probability a of the n estimations 
falls within the correct bin on average, while Pn estimations fall 
outside the correct bin. The standard deviation of the 
distribution may therefore be expressed as 

STD"" (T) = aSTD"" (r^) + /XSTP" (Z") = a^l ^fial. (26) 
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In this analysis, the accuracy breakpoint is defined as the 
noise level for which the contribution of T*^ to the total error 
becomes dominant. It is noted that the RMSE is significantly 
greater than that given by the uniform distribution on. Ia alone 
5 when 



2 . 2 

ari +cr„ 



(27) 



The probability breakpoint may therefore be defined as 

10 

Oj =_i_2 ^^—r-. (28) 



Accordingly, the SNR corresponding to the accuracy breakpoint may 
be determined as the SNR value for which equation (25) equals 

15 equation (28) . 

As indicated above, the method of the second embodiment of 
the sonar system 100 (see Fig. la) includes calculating the mode 
of the echo delay estimation distribution generated from raw or 
preprocessed sonar signals. Such calculation of the distribution 

20 mode allows this second embodiment of the system 100 to increase 
its tolerance to noise. This will be better understood by 
reference to the following analysis. 

A traditional way of combining infoirmation from multiple 
observations is to perform an averaging operation. Because the n 

25 echo delay estimations described above are independent and 
identically distributed, the central limit theorem indicates that 
the standard deviation (i.e., the error) of the averaged random 

vector is times smaller than the error corresponding to each 

of the n separate estimations. This holds true for the region of 
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the peak variability curve (e.g., the peak variability curves 301- 
307) corresponding to the coherent receiver. However, the 
averaging operation does not allow the accuracy breakpoints of the 
curves 301-307 to be shifted to lower SNR values. Accordingly, 
5 while the averaging operation may improve accuracy, it generally 
does not increase noise tolerance. 

For example, consider sampling from a unifoinu distribution 
within an interval Ia with probability a, and sampling from a 
uniform distribution within an interval Iq and probability p, in 
10 which the interval Iq includes a gap corresponding to the interval 
Ia. In the event the distributions are sampled n times to obtain 
Ti,T2f...#Tn samples, an estimate for the time delay may be expressed 
in terms of the sample mean, i.e., 

15 r=-Er,. (29) 

On average, an values ^T^ ,>..yT^ are within the correct bin, 
and pn values yT^ y*^*^T^ are disposed outside of the correct bin. 
The sample mean may therefore be decomposed into two parts, i.e.. 




(30) 



Applying the central limit theorem to the sums in equation (30) 
yields 

25 

.2 



cina\ 



(31) 
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in which '''' STD/' denotes the standard deviation of the delta-like 
distribution (see Fig. 10), and STD^" denotes the standard 
deviation of the uniform distribution (see Fig. 11) - It is noted 
5 that the root mean square error (RMSE) is significantly greater 
than that given by the delta distribution alone when final ^ r 
i.e., 

10 

It is noted that the bound on the probability expressed in 
equation (32) is equal to the bound on the probability expressed 
in equation (27) for a single ping, thereby indicating that the 
bound on the probability a does not change with the niamber of 

15 pings. For this re.ason, taking the mean of the distribution of 
echo delay estimations, as described above with reference to the 
first embodiment of the presently disclosed sonar system, 
generally does not improve the noise tolerance of the system, even 
though the accuracy of echo delay estimation is improved. 

20 The method of the second embodiment of the sonar system ICQ 

(see Fig. la) increases the accuracy of echo delay estimation and 
improves the noise tolerance of the system via a calculation of 
the mode of the distribution of echo delay estimates derived from 
multiple pings. As depicted in Fig. 12, the a priori window of 

25 the cross correlation may be divided into a total of /n = [2X/Aj 

intervals 81,62/-/ Bm, in which Bi=Ia represents the correct bin. 
If PirP2/.../Pm denote the probabilities of estimates falling in each 
of the intervals Bi,B2,...,Bm/ respectively, and Yi, Y2,.-.,Yxa are random 
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5 



variables representing the number of estimates falling in each 
interval, then 



sr. =n, 

(33) 



The joint probability distribution function for the echo delay 
estimates within each bin may be expressed as 



10 

In this second embodiment, the probability of choosing the correct 
bin after calculating the mode of the distribution is defined as 
the probability that the number of estimates disposed within the 
correct bin ki is greater than the number of estimates disposed in 
15 any other bin ki, i^^l, i.e.. 



It is noted that the sum included in equation (35) may be 
20 decomposed into two parts, namely, the probability P^sm^ that more 
than half of the ^^n" estimations fall into the correct bin, and 
the probability jP<5o% that even if less than half of the n 
estimations fall into the correct bin, the number of estimations 
in the correct bin is greater than that of any other bin. 
25 Equation (35) may therefore be expressed as 
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^correct -^>50% '*"-P<50%> (36) 



in which 



15 



20 



p = 



In the second embodiment^ the probability p of an echo delay 
estimate falling outside of the correct bin is substantially 
uniform over the length L of the a priori window. The probability 
10 of an echo delay estimate falling within any given bin or interval 
A of the a priori window is therefore equal to p/(m-l), and 

Pi = a, 

Pj = (l-a)/(m-l), j 9t 1, (38) 



in which ^'Pi" is the probability corresponding to the correct bin 
and '^Pj'' is the probability corresponding to all other bins other 
than the correct bin. Substituting equations (38) into equation 
(37), 

P>smi niay be expressed as 



^>«^ ~t>?/2 



ccy'\ (39) 



in which ^'a'' is a function of the SNR of the environment (see 
equation (25)). The SNR for which P^sm'='<^o (see also equation 
25 (28)) represents an upper bound on the accuracy breakpoint for 
echo delay estimations calculated using the mode of ^^n'' estimates, 
i.e.. 
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a, 




n-k 



(40) 



It is noted that the upper bound corresponding to the total 



probability of choosing the correct bin, as expressed in equation 
5 (36), is lower than that expressed in equation (40), i.e.. 



The accuracy breakpoint that can be achieved via the calculation 

10 of the mode of the echo delay estimation distribution therefore 
corresponds to a lower SNR than that expressed in equation (41), 
which is significantly better than the accuracy breakpoint that is- 
achievable from a single ping, or from a calculation of the mean 
of the distribution of echo delay estimations derived from 

15 multiple pings. 

The second embodiment of the presently disclosed sonar 
system 100 (see Fig. la) will be better understood by reference to 
the following illustrative example. In this example, 

distributions of echo delay estimations are obtained via Monte 

20 Carlo simulations using a cosine packet. Figs. 13a-13d depict 
simulated distributions 1300a-1300d of echo delay estimations in 
histogram form for different SNR values. As shown in Fig. 13a, 
for high SNR values (i.e., ^.20 dB) , essentially all of the echo 
delay estimates are disposed within the central or correct bin. 

25 As the SNR decreases (i.e., as the noise level increases), 
significant errors in the echo delay estimates appear. As shown* 
in Figs. 13b-13d, such errors are substantially uniformly 
distributed over the length of the a priori window, and the ratio 
of the correct'' echo delay estimates (corresponding to the 

30 central bin) and the level of the uniform error distribution 



SNR^n < SNR 



>50%* 



(41) 
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outside of the central bin generally decreases with the level of. 
SNR. It is noted^ however, that even in the presence of 
relatively high levels of noise, the peaks of the distributions 
within the central or correct bins are significantly greater than 
5 those within any one of the other bins. 

Figs. 14a-14d depict representative peak variability curves 
corresponding to different numbers of pings. Specifically/ Fig. 
14a depicts a peak variability curve 1402 for a single ping. As 
shown in Fig. 14a, for high SNR values (e.g., greater than about 

10 17 dB) , the accuracy of echo delay estimation is consistent with 
the performance of a coherent receiver. For lower SNR values 
(e.g., less than about 17 dB) , there is an accuracy breakpoint 
indicating a sharp transition from low RMSE levels to high RMSE 
levels. Figs, 14b-14d depict peak variability curves 1404b-1404d, 

15 1406b-1406d, and 1408b-1408d including accuracy breakpoints for 
multiple pings, e.g., 10 pings (see Fig. 14b), 50 pings (see Fig. 
14c) , and 100 pings (see Fig. 14d) . Specifically, the curves 
1404b-1404d and 1406b-1406d are derived from calculations of the 
mean and the mode, respectively, of the distribution of echo delay 

20 estimates. Further, the curves 1408b-1408d correspond to optimal 
cross correlations of the multiple echo/ping pairs. As shown in 
Figs. 14b-14d, the accuracy breakpoints included in the curves 
140Bb-1408d represent the optimal breakpoints for a stationary 
sonar system and a stationary target. The optimal accuracy 

25 breakpoints are difficult to obtain, however, due to the need for 
precise registration of the multiple echo /ping pairs. Such 
precise registration of echoes/pings normally requires the 
distance between the sonar system and the target to be either 
constant or known in advance, which is not always the case in 

30 practical system applications. 

As indicated above. Figs. 14b-14d also depict the peak 
variability curves 1404b-1404d and 1406b-1406d resulting from the 
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calculations of the mean and the mode, respectively, of the echo 
delay estimation distribution for multiple pings. As shown in 
Figs. 14b-14d, the curves 1404b, 1406b, 1408b, and 1410b, the 
curves 1404c, 1406c, 1408c, and 1410c, and the curves 1404d, 
5 1406d, 1408d, and 1410d gradually shift to lower RMSE values as' 
the number of pings increases, thereby indicating increased 
accuracy. Further, the accuracy breakpoints included in the 
curves 1406b-1406d corresponding to the mode calculation are 
shifted to lower SNR values as the number of pings increases, 

10 thereby indicating increased noise resiliency. No such shifting 
occurs for the accuracy breakpoints included in the curves 1404b- 
1404d corresponding to the calculation of the mean of the echo 
delay estimation distribution. 

Fig. 15 depicts a plurality of performance curves 1504-1508 

15 comprising composites of the peak variability curves 1404b-1404d, 
1406b-1406d, and 1408b-1408d, respectively, showing the accuracy 
breakpoints versus the number of pings. Specifically, the curve 
1508 corresponds to the optimal breakpoints for the cross 
correlation of multiple echo/ping pairs, the curve 1504 

20 corresponds to the accuracy breakpoints obtained after calculating 
the mean of the distribution of echo delay estimates, and the 
curve 1506 corresponds to the accuracy breakpoints obtained after 
calculating the mode of the echo delay estimation distribution. 
As shown in Fig. 15, the accuracy breakpoints obtained after 

25 calculating the mode of the distribution are shifted to lower SNR 
values as the number of pings increases, while the accuracy 
breakpoints obtained after calculating the mean of the 
distribution do not change substantially as the niomber of pings 
increases. 

30 Having described the above illustrative embodiments, other 

alternative embodiments or variations may be made. For example,, 
it was described that temporal differences between successive 
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received signals may be estimated by determining the peak arrival 
time differences (PATDs) of the signals. However, it should be 
understood that temporal differences between signals may be 
estimated using any suitable prominent feature of the signals such 
5 as prominent peaks, valleys, energies, and/or zero crossings of 
the signals. 

It was also described with reference to Fig. la that the 
receiver 104 provides the received signals directly to the time 
delay estimator 108. However, Fig. la merely depicts an 

10 illustrative embodiment of the system 100 and other alternative 
embodiments or variations may be made. For example. Fig. lb 
depicts a system 100b that includes multiple receivers 104. 1-104. n 
providing indications of received signals to a beamformer 105, 
which in turn provides beams to the time delay estimator 108. For 

15 example, the beamformer 105 included in the system 100b may be 
adapted to produce seismic sonar beams in a seismic sonar system. 

It was also described that statistical estimates of time 
delay may be determined by calculating the mean or the mode of the 
distributions of time delay estimations. However, in alternative 

20 embodiments, statistical estimates of time delay may be determined 
by calculating the median of the distribution of time delay 
estimations. Figs. 14b-14d depict illustrative peak variability 
curves 1410b-1410d resulting from calculations of the median of 
the time delay estimation distribution. Like the curves 1406b- 

25 1406d corresponding to the mode calculation, the accuracy 
breakpoints included in the curves 1410b-1410d are shifted to 
lower SNR values as the number of pings increases. However, the 
accuracy breakpoints included in the curves 1406b-1406d are 
shifted to lower SNR values than the breakpoints included in the 

30 curves 1410b-1410d, thereby indicating better noise resiliency for 
the mode calculation. Fig. 15 depicts a performance curve 1510 
comprising a composite of the peak variability curves 1410b-1410d 
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corresponding to the median calculation. As shown in Fig. 15, the 
calculation of the median of the distribution of time delay 
estimates obtained from multiple pings shifts the accuracy 
breakpoints to SNR levels that fall between the SNR levels 
5 obtained from the mean and mode calculations. 

Increased accuracy in time delay estimation was also 
illustrated via improved discrimination of a target comprising a 
cylindrical container filled with different liquids such as fresh 
water, saline water, glycerol, or kerosene. It should be 

10 understood, however, that variations of this illustrative 
technique may be employed. For example, in the field of medical 
ultrasound, multiple pings may be transmitted toward a region of a 
patient's body to determine bone density. When each ping, 
penetrates the bone^ a prominent feature such as the signal peak 

15 is maintained after being bounced from the front and the back of 
the bone. Two peaks can therefore be detected in successive 
echoes, and the time difference between the two peaks can be 
measured. Accordingly^ the patient's bone density can be 
determined by analyzing the temporal differences between 

20 successive echoes. 

Similarly, multiple pings may be transmitted toward the 
patient's heart to determine characteristics of the heart such as 
a heart wall thickness. When each ping penetrates the heart wall, 
a prominent feature such as the signal peak is maintained after 

25 being bounced from the front and the back of the heart wall, and 
the time difference between two signal peaks can be measured. 
Accordingly, the patient's heart wall thickness can be determined 
by analyzing the temporal differences between successive echoes. 
In addition, three dimensional representations of the patient' s 

30 bones, heart, and any other suitable organ may be obtained by 
further application of this illustrative technique. 
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It was also described that the transmitter of the presently 
disclosed embodiment is configured to transmit one or more pings 
through a transmission medium such as water, and that the pings 
travel through the water until they strike an object or target in 
5 the water, which returns one or more echoes toward the sonar 
sensors. Because such targets typically comprise many 

substructures, multiple echoes are typically generated for each 
ping. It should be noted that the size of a window around each 
possible maximum of the envelope of an echo can be determined to 

10 distinguish between peaks representing a density change inside the 
target and - spurious peaks resulting from noise. For example, the 
density change inside the target may correspond to layer 
penetration in geological exploration. When there is mutual 
movement between the transmitter/receiver and the target, such 

15 window determination has to be accomplished for multiple echoes 
from each pinging source separately. In one embodiment, such 
window determination is accomplished as follows. 

First, motion estimation is applied for each set of echoes. 
The motion estimation is based on the assumption that the motion 

20 between the sensor array and the target is rigid, thereby 
utilizing all sensor recordings for each ping to deteantiine the two 
dimensional displacement and rotation of the sensor array relative 
to the target. Because the sensor array is rigid, the same echo 
arriving at all of the sensors forms a noisy straight line in a 

25 plane defined by the sensor axis and an axis corresponding to the 
time of arrival of the echo at each of the sensors. The time of 
arrival or the peak energy falls roughly on a straight line due to 
noise. From the set of points, the properties of the line can be 
estimated, namely, the displacement from the origin and tilt. 

30 This corresponds to the motion estimation of the sonar for each 
ping. A transformation is then performed so that all lines 
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corresponding to echoes from different pings fall on the same 
line. In this way, motion correction is achieved. 

Next, each set of echoes resulting from a single ping is 
shifted to be aligned with echoes resulting from the other pings. 
5 When the properties of the line are estimated, the estimated tilt 
corresponds to the tilt of the sensor array relative to the 
target. Each set of echoes is shifted by shifting the time of 
arrival for each sensor based on a foimula ^'ax+b^', in which ^^x''' 
corresponds to the sensor number 1, 2, etc., ^'a" corresponds to 
10 the tilt ratio, and ^^b" corresponds to a global shift of the 
sensor. 

The analog mean of the energies of all signals is then 
calculated. This constitutes a first estimation of the sonar 
image from the multiple pings. It is noted that the mean signal 

15 is taken with respect to echoes resulting from different pings. 
For each ping, a time series is determined corresponding to the 
echoes. Such a time series is obtained at each sensor. The 
motion estimation and correction described above results in a 
shift in time of the time series, so that for all sensors the time 

20 series will be aligned, and for all pings the time series of each 
sensor will be aligned. Then, for each sensor, an average of 
energy over all pings is determined. 

Next, a threshold is estimated from the analog mean image to 
determine a certain percentile of signal energy defining a window 

25 region around possible energy distribution peaks. Thus, for each 
sensor, there is a single time series corresponding to the average 
over the multiple pings. The standard deviation (STD) of the 
energy of the time series can then be calculated, and the 
threshold can be estimated to be about 3*STD- This threshold may 

30 be employed to define multiple windows around peaks that exceed 
3*STD. The window width corresponds to the region around such 
peaks . 
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Based on the estimated threshold^ temporal regions are then 
defined around each possible energy distribution peak. Finally, 
the median, the mode, or any other suitable statistical estimate 
is calculated from the distribution within each window region. 
5 The separation into such temporal windows constitutes a first 
phase of outliers removal for each peak dete2niiination. 

It will be further appreciated by those of ordinary skill in 
the art that modifications to and variations of the above- 
described improved echo delay estimates from multiple sonar pings 
10 may be made without departing from the inventive concepts 
disclosed herein. Accordingly, the invention should not be viewed 
as limited except as by the scope and spirit of the appended 
claims . 
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